Ab initio mechanical response: internal friction and structure of divacancies in silicon. 
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This letter introduces ab initio study of the full activation- volume tensor of crystalline defects as a 
means to make contact with mechanical response experiments. We present a theoretical framework 
for prediction of the internal friction associated with divacancy defects and give the first ab initio 
value for this quantity in silicon. Finally, making connection with defect alignment studies, we give 
the first unambiguous resolution of the debate surrounding ab initio verification of the ground-state 
structure of the defect. 
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Among the most powerful applications of first prin- 
ciples ab initio approaches in condensed matter sys- 
tems is the interpretation of experimental signatures 
from defects. The extremely efficient, albeit approx- 
imate, functionals available for density-functional the- 
ory have given this approach wide application in con- 
junction with experimental probes such as scanning 
tunneling microscopy, electron paramagnetic resonance, 
electron-nuclear double resonance, and nuclear magnetic 
resonance 0, EI1I1I1E10,0I1GI1III1 Macroscopic 
mechanical response experiments also have proved an ex- 
tremely valuable tool in the study of defects in solid-state 
systems [H 013110, providing key information 
on such issues as defect symmetries and concentrations. 
However, mechanical response studies remain largely ig- 
nored by the ab initio community to date. 

Despite their successes, the aforementioned density- 
functional studies suffer a fundamental flaw: no underly- 
ing theorem ensures that density-functional theory pro- 
vides the energy- or angular-momentum- resolved den- 
sities probed in the experiments involved in the appli- 
cation. Beyond this matter of principle, such quantities 
are not among those which available approximate func- 
tionals predict most reliably. This letter notes that the 
most reliable quantities which density-functional theory 
predicts (bond lengths, bond angles, and lattice param- 
eters) relate directly to the key coupling parameter in 
macroscopic mechanical response experiments, the full 
activation-volume tensor. We thus propose ab initio 
study of mechanical response functions associated with 
this tensor as a powerful and particularly reliable new 
tool in the study of defects in condensed matter systems. 

One such response function, internal friction, is 
a topic of current interest both experimentally and 
theoretically [H H H HO, HJ 0. Below we develop 
the theory of internal friction as applied to divacancy 
defects to provide the first parameter-free ab initio de- 
termination of friction from a point defect, in this case 
the singly negatively charged divacancy in silicon (Si- 
V 2 ~). Although the structure of this defect has been 
inferred through a combination of experimental signa- 
tures and symmetry of electronic states 23J, ab initio 
work to confirm the structure has resulted in an ongoing 



debate[6|, |2J, l'2al |26( which has arisen ultimately from 
a focus on delicate energy differences at or beyond the 
limits of accuracy of current density functionals. Here, 
we demonstrate the power of working with mechanical 
response functions by providing a clear signature which 
resolves the debate. 

Mechanical response of defects — Within linear re- 
sponse, the stress-dependent part of the energy of a point 
defect AE has the form 



AE = 



•J > 



(1) 



where we employ repeated index summation notation 
here and throughout this work and where a and A are, 
respectively, the externally applied stress tensor and the 
defect activation-volume tensor. 

The activation- volume tensor A is accessible in terms of 
equilibrium supercell lattice constants, which are among 
the quantities most reliably determined within density- 
functional theory. This connection comes directly from 
the principle that the strain u induced in a crystal con- 
taining a concentration n of defects with activation vol- 
ume tensor A is = nAy, a result of minimizing the 
sum of the bulk and defect elastic energies. Accordingly, 
one may determine A simply by creating a supercell con- 
taining a single defect and relaxing the defect structure 
along with the supercell lattice vectors. 

Under applied stress, defects with lower symmetry 
than the host crystal tend to reorient so as to minimize 
the energy. Experimental stress-alignment studies, which 
observe the relative thermalized populations of different 
orientations of specific defect types as a function of ap- 
plied stress, thus allow direct access to certain linear com- 
binations of Aij for each type |2^|. Such thermalization 
of defect populations under time-varying external stress 
o~ij(t) provide a mechanism for internal friction, dissipa- 
tion of mechanical energy throughout the bulk of a mate- 
rial. We review this process here in some depth because, 
although our overall logic is the same, the final result 
for the divacancy differs from the oft-quoted result [2?| . 
which assumes that all defect orientations relax among 
each other at equal rates and hence does not apply to 
divacancies. 
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From Q, the total energy per unit volume stored at 
time t among all defect orientations m of a specific defect 
type is 



AE(t) 



(2) 



where n is the total number density of defects of this 
type, and P m and A£™ are the probability and activation 
volume tensor associated with each orientation. Dissipa- 
tion results from energy lost irreversibly to the heat bath 
through transitions among defect orientations and thus 
occurs at the rate 



dE_ 
~dt 



,^ dP n 



A>i 3 -(t). 



(3) 



Finally, stresses ultimately drive the transitions through 
the master equation 

dP m / \ / 

"^T = £ «W (l + /3A% <7«(t)J P™ (t) 
m' 

-i/ w (l + /3A7; ( x y -(t))P m (t), (4) 

where we take transitions from orientation m' to m to 
be thermally activated with rate v mm i in the absence of 
external stress and with stress dependencies which we 
have linearized. Solving an d substituting the result 
into ® then gives the final dissipation rate. 

Under the not uncommon special case underlying the 
result |27j , in which the zero-stress transition rates among 
different orientations of a defect type are equal by sym- 
metry, v mm ' = v, the master equation may be solved 
analytically, resulting in a defect contribution to the in- 
verse quality factor Q _1 , the fraction of energy lost per 
radian of oscillation phase, of 



Q- 1 = n (3- 
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(5) 
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Re {dij exp (iwt)), AAJ 
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g Em' -^-ii ' i s the standard elastic compliance four- 
tensor, and Lij-ki is an anelastic four-tensor sharing the 
symmetry of the full set of defect orientations. 

In a system with more than one type of defect, the 
inverse quality factors for each type add. If the transition 
rate v a associated with each type s has the same value v, 
then Q -1 can be written as JHJ) with 



j tot 



-L 



ij ; kl 



(6) 



replacing Lij.^i, where n s is the number density for each 
type of defect so that ^ s n s = n. Below we show how to 
apply Eqs. lloKit to find the correct result for divacancies. 



Finally, we note here that there is a simple, exact rela- 
tionship between the anelastic four-tensor and the results 
from experimental studies of the inverse quality factor 
of a mechanical oscillator as a function of temperature, 
which generally show a peak when the thermally acti- 
vated, and thus temperature-dependent, transition rate 
v(T) corresponds to the oscillator frequency [l2j. In par- 
ticular, from (!')!( ill we have that the experimentally ac- 
cessible quantity fceTQ -1 , where fc# is Boltzmann's con- 
stant, has a maximum at precisely v(T) = uj / g with value 



max (ksTQ 1 ) = 



~ * T tOt ~ 

2 V*jSij;klVkl 
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Ab initio determination of L*°* fe; thus gives a parameter- 
free relationship between the maximum temperature- 
internal friction product and the defect concentration, 
thereby allowing mechanical response experiments to 
provide a direct, parameter- free measure of absolute de- 
fect concentrations for the first time. 

Application to divacancy in silicon — The importance 
of the divacancy in silicon has led to its extensive study, 
both theoretically and experimentallyH IH |H Hi 
113,01. In contrast to single vacancies, which are quite 
mobile and thus anneal readily, divacancies have low mo- 
bility and are among the most common stable defects in 
silicon at room temperature. The defect has four charge 
states, singly positively charged, neutral, singly and dou- 
bly negatively charged, with the singly charged defect 
(Si-KT) playing an important role in carrier recombina- 
tion o. 

Since the publication of Watkins and Corbett's pio- 
neering study (2^, the ab initio determination of the 
precise nature of the ground-state structure of Si-V^ - has 



become the subject of debate 6, 24, 25, 26] . The idealized 
defect has D^d symmetry along the (111) axis connect- 
ing the sites of the neighboring vacant atoms. The de- 
fect also introduces partially filled degenerate electronic 
states into the gap and thus undergoes a Jahn- Teller dis- 
tortion which ultimately lowers the symmetry to C2/1. 
The debate arises because two stable structures, termed 
pairing and resonant, with very similar energies are con- 
sistent with the C^h symmetry of the defect. 

Figure n shows the two competing ground-state struc- 
tures as projected along the (111) defect axis connecting 
the sites of the two vacant atoms. The pairing configura- 
tion breaks symmetry by moving two pairs of atoms (ah 
and a'b' in the figure) toward each other along a (110) 
reconstruction axis to form stronger bonds at the expense 
of strain energy. In the resonant structure, the same pairs 
of atoms move away from each other to form a less favor- 
able bonding configuration at a nearly correspondingly 
lower cost in strain. Whenever this work requires a spe- 
cific coordinate system, it shall be such that the defect 
axis is [111] and the bonding axis is [110]. 
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FIG. 1: Projection of divacancy structure along (111) sym- 
metry axis: missing atoms (dashed circles), atoms in plane 
closer to observer (larger solid circles), atoms in plane further 
from observer (smaller solid circles). 

The literature reports a number of density-functional 
values for the energy difference AE = E pa i r — E les be- 
tween the two competing configurations: 0.0024 eV[24j. 
~ eV|H, -0.12 eVfc], -0.18 eV (present work). The 
fact that these differences are all quite small and of the 
order of the uncertainties due to the approximate density 
functionals which these works employ (all other compu- 
tational uncertainties notwithstanding!) underscores the 
difficulty of using total energies to resolve the ground 
state structure of the Si-V 2 _ defect and indicates that 
previous studies are inconclusive. We propose instead to 
use the considerably different activation volume tensors 
of the competing reconstructions as a more appropriate 
signature to confirm the ground state. 

Experimentally, Watkins and Corbett |2^| explored the 
activation volume tensor of the Si-V 2 _ defect in depth us- 
ing electron spin resonance to study thermal alignment 
of defect subpopulations under external (110) stresses. 
These experiments, as do also internal friction experi- 
ments, take place under conditions where (110) recon- 
struction axes have time to thermalize while (111) defect 
axes do not. Such (110) stresses split divacancies into 
two classes, a and /3, according to the orientation of the 
defect axis, with a corresponding to the defect axis be- 
ing perpendicular to the stress. Within each of these 
classes, there is a further splitting of the energy into two 
distinct values, for a total of four energetically different 
states [33). 

As thermalization occurs only among choices of recon- 
struction axis, the quantities which stress-alignment ex- 
periments actually access are the energy splittings within 
each class, AE a and AEp, respectively, each of which re- 
late directly to certain linear combinations of components 
the defect activation- volume tensor, 

AE a = ~ (A n + 2A 12 - 2A 13 - A 33 ) = -ctAA q (8) 
AEp = -| (A n - 2A 12 + 2A 13 - A 33 ) = -ctAA^, 

where a is the magnitude of the external (110) stress and 



Aij are the Cartesian components of the activation vol- 
ume tensor in the cubic coordinate system defined above. 

Calculations and results - The ab initio elec- 
tronic structure calculations below employ the total- 
energy plane-wave density-functional pseudopotential 
approach |34j | within the local spin-density approxima- 
tion (LSDA) using a pseudopotential of the Kleinman- 
Bylander form J35] with p and d non-local corrections. 
The calculations expand the Kohn-Sham orbitals in a 
plane-wave basis set with a cutoff energy of 6 Hartrees 
within a cubic sixty-four atom supercell, sampling the 
Brillioun zone at eight k-points, reduced to four by 
time reversal symmetry. Finally, we employ the analyti- 
cally continued functional approach pjif to minimize the 
Kohn-Sham energy with respect to the electronic degrees 
of freedom. 

To determine A, we relax a supercell containing a sin- 
gle defect and compute the relaxed strain and thus A. 
In principle, this can be done in a single joint relaxation 
with the internal atomic coordinates. For the present 
study, however, we followed the nearly equivalent ap- 
proach of minimizing the internal coordinates while mini- 
mizing separately along each of the six independent com- 
ponents of supercell strain while including the Poisson 
ratios of the supercell when appropriate. 

Ogiit and Chelikowskyj^ emphasize the need to tailor 
supercell shape to accommodate the relaxation pattern 
of V% to obtain accurate results for defect energy differ- 
ences. To explore supercell-size effects on the extraction 
of the activation-volume tensor, we carried out a con- 
vergence study of AA a and AA^ from JHJ employing the 
environment-dependent interatomic potential (EDIP) for 
silicon ^37] for cubic supercells of lattice constant from 
2ao through 6ao, where oq is the lattice constant of the 
"primitive" eight-atom cubic cell. Over this range of cell 
sizes, we observe a total change of only 12.5% (6.3%) in 
AA Q (AA^), with 70% (86%) of the change occurring 
between 2ao and 3oq. We thus conclude that A is not a 
particularly sensitive function of cell size and that a su- 
percell of 64 atoms suffices to give ab initio values with 
an uncertainty on the order of 10%. 

Table [I] summarizes our ab initio results for the activa- 
tion volume tensors of the two candidate defect struc- 
tures, and Figure El compares our predictions directly 
with the experimentally available linear combinations, 
AA Q and AA^ [2^. Our predictions for the resonant 
configuration are clearly inconsistent with the measure- 
ments, whereas our results for the pairing configuration 
show good agreement with errors (+20% and -6% for 
AA Q and AA^, respectively) consistent with supercell- 
convergence uncertainty we estimate above. The figure 
also compares our ab initio prediction of another linear 
combination of components of the activation volume ten- 
sor, B33 = C33 ; ijAij where Cij-^i is the elastic constant 
four-tensor, with an estimate of this quantity from |23l |. 

Turning finally to the internal friction, because the de- 
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-11 

18 
1.5 



18 1.5 
-11 1.5 
1.5 10 



0.6 -1.5 
-1.5 0.6 
10 10 



10 
10 

-13 



pair res 



Ln 100 38 
L12 -50 -19 
L 44 2 9 3 2 79 



TABLE I: Activation-volume and anelastic tensors for com- 
peting ground-state structures of singly negatively charged 
divacancy in silicon 
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FIG. 2: DFT results and comparison to experiment of AA a , 
AA/3 and B33. For each data set the black bar corresponds 
to the experimental data, the white bar to our DFT results 
for the pairing configuration and the gray to those for the 
resonant configuration. B33 has been scaled by Cu for display 
purposes. 



feet axes do not thermalize under typical experimental 
conditions, each choice of defect axis must be treated as 
a separate type of defect s in (JTjJ). Under normal sam- 
ple preparation, all four distinct choices of (111) defect 
axes will occur with equal probabilities, n s — n/4, re- 
sulting in an L tot with cubic symmetry. The thermaliz- 
ing reconstruction axes then constitute the orientations 
m within each type and have sufficient symmetry to en- 
sure the result ©. Table d gives the three unique val- 
ues of the resulting cubic anelastic four-tensor in con- 
tracted notation. This is the first ab initio prediction 
for the components of the anelastic tensor for a defect, 
a quantity of current research interest, particularly in 
the optimization of micro- and nano- electromechanical 
devices. (See, for example 0,H3-) Our value of L tot in- 
dicates that torsional (pure shear) oscillators will experi- 
ence a maximum in divacancy-mediated loss in the range 
5.5 x l(r 28 cm 3 K < max(TQ- 1 )/n < 8.5 x lCT 28 cm 3 K, 
with the precise value depending on the crystallographic 
orientation of the device through J^J. 

In conclusion, this letter introduces ab initio study 
of the full activation-volume tensor of crystalline de- 
fects as a quantity which current approximate density- 
functionals give accurately and which is of direct use in 
making contact with mechanical response experiments, 
including both stress-alignment studies and measure- 
ments of macroscopic internal friction. Illustrating the 



power of the approach, this letter gives the first unam- 
biguous ab initio verification of the nature of the ground 
state of the singly negatively charged divacancy in sili- 
con and the first parameter-free theoretical calculation of 
the peak internal friction associated with a point defect. 
This latter quantity then forms the basis for a straightfor- 
ward method for determining defect concentrations via ab 
initio interpretation of macroscopic mechanical response 
experiments. 

Primary support of this work is through the MRSEC 
program of the NSF (No. DMR-0079992). 
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